# Load Packages
library(frontier)
library(quantreg)
library(npbr)
library(Rglpk)
library(data.table)
library(RColorBrewer)
library(zoo)
library(pracma)
library(plot3D)
library(ggplot2)
library(memsic)
library(xtable)
library(countrycode)

boundary_routine <- function(dat, input = input, output = output, method)
{
  dat$x <- input
  dat$y <- output
  id <- order(dat$x)
  divide_by_this <- (max(dat$x) - min(dat$x)) * (max(dat$y) - min(dat$y))
  ymin <- min(dat$y)
  idx <- sample(1:nrow(dat), nrow(dat), replace=TRUE)
  sample.dta <- dat[idx,]
  y <- sample.dta$y
  x <- sample.dta$x
  if(method == "Polynomial")
  {
    odeg <- opt_degree(x, y, x, prange=0:20)
    polfront <- poly_estimate(x, y, x, deg = odeg)
    AUC_poly <- sum(diff(sort(x)) * rollmean(sort(polfront), 2))
    AUC_lower <- sum(diff(dat$x[id]) * ymin)
    AUC_difference <- AUC_poly - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    result <- 1 - AUC_percent
  }
  if(method == "Kernel")
  {
    bw <- kern_smooth_bw(x, y)
    kernsmooth <- kernel_smoothing(x, y, x, h = bw)
    AUC_kern <- sum(diff(sort(x)) * rollmean(sort(kernsmooth), 2))
    AUC_lower <- sum(diff(dat$x[id]) * ymin)
    AUC_difference <- AUC_kern - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    result <- 1 - AUC_percent
  }
  if(method == "SFA")
  {
    ka_sfa <- sfa_est(y ~ x | x, ineffDecrease = TRUE, data = sample.dta)
    y_sfa <- ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * x
    AUC_sfa <- sum(diff(sort(dat$x)) * rollmean(sort(y_sfa), 2))
    AUC_lower <- sum(diff(dat$x[id]) * ymin)
    AUC_difference <- AUC_sfa - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    result <- 1 - AUC_percent
  }
  if(method == "QR")
  {
    ka_qr <- rq(y ~ x, tau = 0.95, data = sample.dta)
    y_qr <- ka_qr$coefficients[1] + ka_qr$coefficients[2] * input
    AUC_qr <- sum(diff(sort(sample.dta$x)) * rollmean(sort(y_qr), 2))
    AUC_lower <- sum(diff(sort(sample.dta$x)) * ymin)
    AUC_difference <- AUC_qr - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    result <- 1 - AUC_percent
  }
  result
}

which_technique_to_use <- function(method, dat, input = input, output = output)
{
  id <- order(input)
  divide_by_this <- (max(input) - min(input)) * (max(output) - min(output))
  ymin <- min(dat$y)
  if(method == "QR"){
    ka_qr <- rq(output ~ input, tau = 0.95, data = dat)
    y_qr <- ka_qr$coefficients[1] + ka_qr$coefficients[2] * input
    AUC_qr <- sum(diff(input[id]) * rollmean(sort(y_qr), 2))
    AUC_lower <- sum(diff(input[id]) * ymin)
    AUC_difference <- AUC_qr - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    qr_result <- 1 - AUC_percent
    return(qr_result)
  }
  if(method == "SFA"){
    ka_sfa <- sfa_est(output ~ input | input, ineffDecrease = TRUE, data = dat)
    y_sfa <- ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * input
    AUC_sfa <- sum(diff(input[id]) * rollmean(sort(y_sfa), 2))
    AUC_lower <- sum(diff(input[id]) * ymin)
    AUC_difference <- AUC_sfa - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    sfa_result <- 1 - AUC_percent
    return(sfa_result)
  }
  
  if(method == "Polynomial"){
    odeg <- opt_degree(xtab = input, ytab = output, x = input, prange = 1:20)
    polfront <- poly_estimate(input, output, x = input, deg = odeg)
    AUC_poly <- sum(diff(input[id]) * rollmean(sort(polfront), 2))
    AUC_lower <- sum(diff(input[id]) * ymin)
    AUC_difference <- AUC_poly - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    poly_result <- 1 - AUC_percent
    return(poly_result)
  }
  
  if(method == "Kernel"){
    bw <- kern_smooth_bw(input, output, method = "u", technique = "noh",
                         bw_method="bic")
    kernsmooth <- kernel_smoothing(input, output, input, h = bw)
    AUC_kern <- sum(diff(input[id]) * rollmean(sort(kernsmooth), 2))
    AUC_lower <- sum(diff(input[id]) * ymin)
    AUC_difference <- AUC_kern - AUC_lower
    AUC_percent <- AUC_difference/divide_by_this
    kern_result <- 1 - AUC_percent
    return(kern_result)
  }
}

replicates <- function(method, nboots, dat, input, output)
{
  samples <- replicate(nboots, boundary_routine(dat, method = method,
    input = input, output = output))
}



asciify <- function(df, pad = 1, ...) {
  ## error checking
  stopifnot(is.data.frame(df))
  ## internal functions
  SepLine <- function(n, pad = 1) {
    tmp <- lapply(n, function(x, pad) paste(rep("-", x + (2* pad)),
                                            collapse = ""),
                  pad = pad)
    paste0("+", paste(tmp, collapse = "+"), "+")
  }
  Row <- function(x, n, pad = 1) {
    foo <- function(i, x, n) {
      fmt <- paste0("%", n[i], "s")
      sprintf(fmt, as.character(x[i]))
    }
    rowc <- sapply(seq_along(x), foo, x = x, n = n)
    paste0("|", paste(paste0(rep(" ", pad), rowc, rep(" ", pad)),
                      collapse = "|"),
           "|")
  }
  ## convert everything to characters
  df <- as.matrix(df)
  ## nchar in data
  mdf <- apply(df, 2, function(x) max(nchar(x)))
  ## nchar in names
  cnames <- nchar(colnames(df))
  ## max nchar of name+data per elements
  M <- pmax(mdf, cnames)
  ## write the header
  sep <- SepLine(M, pad = pad)
  writeLines(sep)
  writeLines(Row(colnames(df), M, pad = pad))
  writeLines(sep)
  ## write the rows
  for(i in seq_len(nrow(df))) {
    ## write a row
    writeLines(Row(df[i,], M, pad = pad))
    ## write separator
    writeLines(sep)
  }
  invisible(df)
}


desMthing <- function(i, xtab = xtab)
{
  xtab^i
}

gridMthing <- function(i, x = x, xtab = xtab)
{
  x^i
}
poly_estimate <- function(xtab, ytab, x, deg)
{
  opt_coef <- (1/(1:(deg + 1))) * (max(xtab)^(1:(deg + 1)) -
                                     min(xtab)^(1:(deg + 1)))
  desM <- sapply(0:deg, desMthing, xtab = xtab)
  dir <- c(rep(">=", length(xtab)))
  bounds <- list(lower = list(ind = 1:(deg + 1), val = rep(-Inf,
                                                           deg + 1)), upper = list(ind = 1:(deg + 1), val = rep(Inf, deg + 1)))
  Sol <- Rglpk_solve_LP(obj = opt_coef, mat = desM,
                        dir = c(rep(">=", length(xtab))), rhs = ytab,
                        bounds, types = NULL, max = FALSE)
  OPT <- Sol$solution
  gridM <- sapply(0:deg, gridMthing, xtab = xtab, x = x)
  fitt <- gridM %*% OPT
  return(fitt)
}

difs <- function(i, xs, ys, x = x)
{
  log(sum(abs(ys - poly_estimate(xtab = xs,
                                 ytab = ys, x = x, i)))) + (i + 1)/length(xs)
}

opt_degree <- function(xtab, ytab, x, prange = 0:20)
{
  criteria <- cbind(prange, sapply(prange, difs, xs = xtab, ys = ytab, x = x))
  return(prange[which.min(criteria[, 2])])
}

bootstrap_degrees <- function(dta)
{
  opt_degree(xtab = dta$gdp, ytab = dta$life)
}

bs.sample <- function(dta, ...)
{
  idx <- sample(1:nrow(dta), nrow(dta), replace=TRUE)
  dta <- dta[idx,]
  rownames(dta) <- NULL
  dta
}

bootstrap_estimates <- function(i)
{
  dta <- bootstrapped_data[[i]]
  deg <- degrees[[i]]
  poly_estimate(dta[, gdp], dta[, life], x = mydat$gdp, deg = deg)
}

bs.routine <- function(dta)
{
  idx <- sample(1:nrow(dta), nrow(dta), replace=TRUE)
  sample.dta <- dta[idx,]
  y <- sample.dta$y
  x <- sample.dta$x
  odeg <- opt_degree(x, y, analysis_dat$x, prange=0:20)
  polfront <- poly_estimate(x, y, analysis_dat$x, deg = odeg)
  c(1 - polfront)
}

bs.routine.beta <- function(dta)
{
  idx <- sample(1:nrow(dta), nrow(dta), replace=TRUE)
  sample.dta <- dta[idx,]
  y <- sample.dta$y
  x <- sample.dta$x
  bw <- kern_smooth_bw(x, y, method="u", technique="noh", bw_method="bic")
  kernsmooth <- kernel_smoothing(x, y, analysis_dat$x, h = bw)
  c(1 - kernsmooth)
}

bs.routine_general <- function(xtab, dta, method)
{
  idx <- sample(1:nrow(dta), nrow(dta), replace=TRUE)
  sample.dta <- dta[idx,]
  y <- sample.dta$y
  x <- sample.dta$x
  if(method == "Kernel")
  {
    bw <- kern_smooth_bw(x, y, method="u", technique="noh", bw_method="bic")
    kernsmooth <- kernel_smoothing(x, y, x, h = bw)
    return(kernsmooth)
  }
  if(method == "Poly")
  {
    odeg <- opt_degree(x, y, xtab, prange=0:20)
    polfront <- poly_estimate(x, y, xtab, deg = odeg)
    return(polfront)
  }
  if(method == "QR")
  {
    qr <- rq(y ~ x, tau = 0.95, data = sample.dta)
    return(qr$fitted.values)
  }
}

comb_dat <- function(dat, ...)
{
  dat[, gdp]
}

put_wealth_health_together <- function(i)
{
  cbind(wealth[[i]], hw_bs_est[[i]])
}

kern_smooth_bw<-function(xtab, ytab, method="u", technique="noh", bw_method="bic")
{
  n<-length(xtab)
  ndata<-length(xtab) # number of data points
  
  # sorting step to use the Priestly-Chao estimator
  oind<-order(xtab)
  xtab<-xtab[oind]
  ytab<-ytab[oind]
  
  if (technique=="noh" && bw_method=="bic")
  {
    h_min<-2*max(diff(sort(xtab)))
    h_max<-(max(xtab)-min(xtab))
    h_grid<-seq(h_min,h_max,length=184)
    
    BIC<-NULL
    
    for ( h in h_grid)
    {
      Phi<-kernel_smoothing(xtab,ytab,xtab,h,method,technique="noh")
      
      # calculating model complexity
      xtab2<-xtab[-1]
      ytab2<-ytab[-1]
      DX<-diff(xtab)
      
      A<-dnorm(outer(xtab,xtab2,'-')/h)/h
      B<-rep(1,n) %*% t(ytab2*DX)
      S<-A * rep(1,n) %*% t(DX)
      SS<-S[-1,]
      DF<-sum(diag(SS))
      
      BIC<-c(BIC,log(sum(abs(ytab-Phi)))+log(length(ytab))*DF/(2*length(ytab)))
    }
    
    if (which.min(BIC)==1) {mind<-which.min(BIC[-1])+1}
    if (which.min(BIC)>1) {mind<-which.min(BIC)}
    hopt<- h_grid[mind]
  }
  
  if (technique=="noh" && bw_method=="cv")
  {
    bw<-npregbw(xdat=xtab, ydat=ytab, regtype="lc", ckertype="gaussian")
    hopt<-max(max(diff(sort(xtab))),bw$bw)
  }
  
  
  if (technique=="pr" && bw_method=="bic")
  {
    h_min<-2*max(diff(sort(xtab)))
    h_max<-(max(xtab)-min(xtab))
    h_grid<-seq(h_min,h_max,length=30)
    
    BIC<-NULL
    
    for ( h in h_grid)
    {
      Phi<-kernel_smoothing(xtab,ytab,xtab,h,method,technique="pr")
      
      # calculating model complexity
      xtab2<-xtab[-1]
      ytab2<-ytab[-1]
      DX<-diff(xtab)
      
      A<-dnorm(outer(xtab,xtab2,'-')/h)/h
      B<-rep(1,n) %*% t(ytab2*DX)
      S<-A * rep(1,n) %*% t(DX)
      SS<-S[-1,]
      DF<-sum(diag(SS))
      
      BIC<-c(BIC,log(sum(abs(ytab-Phi)))+log(length(ytab))*DF/(2*length(ytab)))
    }
    if (which.min(BIC)==1) {mind<-which.min(BIC[-1])+1}
    if (which.min(BIC)>1) {mind<-which.min(BIC)}
    hopt<- h_grid[mind]
    
  }
  
  if (technique=="pr" && bw_method=="cv")
  {
    bw<-npregbw(xdat=xtab, ydat=ytab, regtype="lc", ckertype="gaussian")
    hopt<-max(max(diff(sort(xtab))),bw$bw)
    
  }
  return(hopt)
  
}

kernel_smoothing <- function (xtab, ytab, x, h, method = "u", technique = "noh")
{
  stopifnot(method %in% c("u", "m", "mc"))
  stopifnot(technique %in% c("pr", "noh"))
  consX <- unique(sort(c(xtab, seq(min(xtab), max(xtab), length = 201))))
  Cstable = 10000
  n <- length(xtab)
  ndata <- length(xtab)
  ncons <- length(consX)
  oind <- order(xtab)
  xtab <- xtab[oind]
  ytab <- ytab[oind]
  if (technique == "noh") {
    oind <- order(xtab)
    xtab <- xtab[oind]
    ytab <- ytab[oind]
    xtab2 <- xtab[-1]
    ytab2 <- ytab[-1]
    DX <- diff(xtab)
    r_end <- max(xtab)
    l_end <- min(xtab)
    opt_coef <- (pnorm((r_end - xtab2)/h) - pnorm((l_end -
                                                     xtab2)/h)) * DX * ytab2
    C0 <- .fitMat(xtab, ytab, xtab, h, type = 0)
    C1 <- .fitMat(xtab, ytab, consX, h, type = 1)
    C2 <- .fitMat(xtab, ytab, consX, h, type = 2)
    obj <- c(opt_coef, -opt_coef)
    if (method == "u") {
      mat_temp <- C0
      mat <- rbind(cbind(mat_temp, -mat_temp), rep(1, length(obj)))
      rhs <- c(ytab, Cstable)
      dir <- c(rep(">=", ndata), "<=")
    }
    if (method == "m") {
      mat_temp <- rbind(C0, C1)
      mat <- rbind(cbind(mat_temp, -mat_temp), rep(1, length(obj)))
      rhs <- c(ytab, rep(0, ncons), Cstable)
      dir <- c(rep(">=", ndata + ncons), "<=")
    }
    if (method == "mc") {
      mat_temp <- rbind(C0, C1, C2)
      mat <- rbind(cbind(mat_temp, -mat_temp), rep(1, length(obj)))
      rhs <- c(ytab, rep(0, 2 * ncons), Cstable)
      dir <- c(rep(">=", ndata + ncons), rep("<=", ncons),
               "<=")
    }
    Sol <- Rglpk_solve_LP(obj, mat, dir, rhs, types = NULL,
                          max = FALSE)
    OPT_temp <- Sol$solution
    OPT <- OPT_temp[1:length(opt_coef)] - OPT_temp[(length(opt_coef) +
                                                      1):(length(opt_coef) * 2)]
    fitt <- .fitMat(xtab, ytab, x, h, 0) %*% OPT
  }
  if (technique == "pr") {
    A <- .fitMat_nw(xtab, ytab, xtab, h, type = 0)
    A.deriv.1 <- .fitMat_nw(xtab, ytab, consX, h, type = 1)
    A.deriv.2 <- .fitMat_nw(xtab, ytab, consX, h, type = 2)
    if (method == "u") {
      p.hat <- solve.QP(Dmat = diag(n), dvec = rep(1, n),
                        Amat = t(A), bvec = ytab, meq = 0)$solution
    }
    if (method == "m") {
      p.hat <- solve.QP(Dmat = diag(n), dvec = rep(1, n),
                        Amat = cbind(t(A), t(A.deriv.1)), bvec = c(ytab,
                                                                   rep(0, ncons)), meq = 0)$solution
    }
    if (method == "mc") {
      p.hat <- solve.QP(Dmat = diag(n), dvec = rep(1, n),
                        Amat = cbind(t(A), t(A.deriv.1), -t(A.deriv.2)),
                        bvec = c(ytab, rep(0, ncons), rep(0, ncons)),
                        meq = 0)$solution
    }
    fitt <- .fitMat_nw(xtab, ytab, x, h, type = 0) %*% p.hat
  }
  return(as.vector(fitt))
}

solve.QP <- function (Dmat, dvec, Amat, bvec, meq = 0, factorized = FALSE)
{
  n <- nrow(Dmat)
  q <- ncol(Amat)
  if (missing(bvec))
    bvec <- rep(0, q)
  if (n != ncol(Dmat))
    stop("Dmat is not symmetric!")
  if (n != length(dvec))
    stop("Dmat and dvec are incompatible!")
  if (n != nrow(Amat))
    stop("Amat and dvec are incompatible!")
  if (q != length(bvec))
    stop("Amat and bvec are incompatible!")
  if ((meq > q) || (meq < 0))
    stop("Value of meq is invalid!")
  iact <- rep(0, q)
  nact <- 0
  r <- min(n, q)
  sol <- rep(0, n)
  lagr <- rep(0, q)
  crval <- 0
  work <- rep(0, 2 * n + r * (r + 5)/2 + 2 * q + 1)
  iter <- rep(0, 2)
  res1 <- .Fortran(.QP_qpgen2, as.double(Dmat), dvec = as.double(dvec),
                   as.integer(n), as.integer(n), sol = as.double(sol), lagr = as.double(lagr),
                   crval = as.double(crval), as.double(Amat), as.double(bvec),
                   as.integer(n), as.integer(q), as.integer(meq), iact = as.integer(iact),
                   nact = as.integer(nact), iter = as.integer(iter), work = as.double(work),
                   ierr = as.integer(factorized))
  if (res1$ierr == 1)
    stop("constraints are inconsistent, no solution!")
  else if (res1$ierr == 2)
    stop("matrix D in quadratic function is not positive definite!")
  list(solution = res1$sol, value = res1$crval, unconstrained.solution = res1$dvec,
       iterations = res1$iter, Lagrangian = res1$lagr, iact = res1$iact[1:res1$nact])
}

.fitMat <- function (X, Y, consX, h, type = 0)
{
  n <- length(X)
  stopifnot(sum((order(X) - 1:n)^2) == 0)
  X2 <- X[-1]
  Y2 <- Y[-1]
  DX <- diff(X)
  A <- dnorm(outer(consX, X2, "-")/h)/h
  B <- rep(1, length(consX)) %*% t(Y2 * DX)
  if (type == 0) {
    return(A * B)
  }
  if (type == 1) {
    A2 <- .dderi(outer(consX, X2, "-")/h)/h^2
    return(A2 * B)
  }
  if (type == 2) {
    A3 <- .d2deri(outer(consX, X2, "-")/h)/h^3
    return(A3 * B)
  }
}

.dderi <- function (x)
{
  return(-x * dnorm(x))
}

.d2deri <- function (x)
{
  return((x^2 - 1) * dnorm(x))
}

sfa_est <- function (formula, data = sys.frame(sys.parent()), ineffDecrease = TRUE, 
                     truncNorm = FALSE, timeEffect = FALSE, startVal = NULL, tol = 1e-05, 
                     maxit = 1000, muBound = 2, bignum = 1e+16, searchStep = 1e-05, 
                     searchTol = 0.001, searchScale = NA, gridSize = 0.1, gridDouble = TRUE, 
                     restartMax = 10, restartFactor = 0.999, printIter = 0) 
{
  formula <- Formula::as.Formula(formula)
  if (length(formula)[2] == 1) {
    modelType <- 1
    effFormula <- NULL
  }
  else if (length(formula)[2] == 2) {
    modelType <- 2
    effFormula <- formula(formula, lhs = 0, rhs = 2)
  }
  else {
    stop("argument 'formula' has an inappropriate number of RHS parts")
  }
  formula <- formula(formula, lhs = 1, rhs = 1)
  if (class(formula) != "formula") {
    stop("argument 'formula' must be a formula")
  }
  else if (length(formula) != 3) {
    stop("argument 'formula' must be a 2-sided formula")
  }
  if (!is.logical(ineffDecrease) || length(ineffDecrease) != 
      1) {
    stop("argument 'ineffDecrease' must be a single logical value")
  }
  if (!is.logical(truncNorm)) {
    stop("argument 'truncNorm' must be logical")
  }
  if (truncNorm && modelType == 2) {
    warning("argument 'truncNorm' is ignored in", " Efficiency Effects Frontiers (EEF)")
  }
  if (!is.logical(timeEffect)) {
    stop("argument 'timeEffect' must be logical")
  }
  if (timeEffect && !"plm.dim" %in% class(data)) {
    warning("argument 'timeEffect' is ignored in case of", 
            " cross-sectional data")
  }
  if (!is.numeric(printIter)) {
    stop("argument 'printIter' must be numeric")
  }
  else if (printIter != round(printIter)) {
    stop("argument 'printIter' must be an iteger")
  }
  else if (printIter < 0) {
    stop("argument 'printIter' must be non-negative")
  }
  printIter <- as.integer(printIter)
  if (length(searchScale) != 1) {
    stop("argument 'searchScale' must be a single logical value or NA")
  }
  else if (is.na(searchScale)) {
    indic <- as.integer(1)
  }
  else if (is.logical(searchScale)) {
    indic <- as.integer(2 - 2 * searchScale)
  }
  else {
    stop("argument 'searchScale' must be a logical value or NA")
  }
  if (!is.numeric(tol)) {
    stop("argument 'tol' must be numeric")
  }
  else if (tol < 0) {
    stop("argument 'tol' must be non-negative")
  }
  if (!is.numeric(searchTol)) {
    stop("argument 'searchTol' must be numeric")
  }
  else if (searchTol < 0) {
    stop("argument 'searchTol' must be non-negative")
  }
  if (!is.numeric(muBound) || length(muBound) != 1) {
    stop("argument 'muBound' must be a numeric scalar")
  }
  else if (is.infinite(muBound)) {
    muBound <- 0
  }
  if (!is.numeric(bignum)) {
    stop("argument 'bignum' must be numeric")
  }
  else if (bignum <= 0) {
    stop("argument 'bignum' must be positive")
  }
  if (!is.numeric(searchStep)) {
    stop("argument 'searchStep' must be numeric")
  }
  else if (searchStep <= 0) {
    stop("argument 'searchStep' must be positive")
  }
  if (!is.logical(gridDouble) || length(gridDouble) != 1) {
    stop("argument 'gridDouble' must be a single logical value")
  }
  if (!is.numeric(gridSize)) {
    stop("argument 'gridSize' must be numeric")
  }
  else if (gridSize <= 0) {
    stop("argument 'gridSize' must be positive")
  }
  if (!is.numeric(maxit) || length(maxit) != 1) {
    stop("argument 'maxit' must be a single numeric scalar")
  }
  else if (maxit != round(maxit)) {
    stop("argument 'maxit' must be an integer")
  }
  else if (maxit < 0) {
    stop("argument 'maxit' must not be negative")
  }
  maxit <- as.integer(maxit)
  if (!is.numeric(restartMax) || length(restartMax) != 1) {
    stop("argument 'restartMax' must be a single numeric scalar")
  }
  else if (restartMax != round(restartMax)) {
    stop("argument 'restartMax' must be an integer")
  }
  else if (restartMax < 0) {
    stop("argument 'restartMax' must not be negative")
  }
  restartMax <- as.integer(restartMax)
  if (!is.numeric(restartFactor) || length(restartFactor) != 
      1) {
    stop("argument 'restartFactor' must be a numeric scalar")
  }
  else if (is.infinite(restartFactor)) {
    stop("argument 'restartFactor' must be finite")
  }
  mc <- match.call(expand.dots = FALSE)
  m <- match("data", names(mc), 0)
  mf <- mc[c(1, m)]
  mf$formula <- formula
  attributes(mf$formula) <- NULL
  mf$na.action <- na.pass
  mf[[1]] <- as.name("model.frame")
  mf <- eval(mf, parent.frame())
  mt <- attr(mf, "terms")
  xMat <- model.matrix(mt, mf)
  xNames <- colnames(xMat)
  yVec <- model.response(mf)
  yName <- as.character(formula)[2]
  if (length(yVec) != nrow(xMat)) {
    stop("the number of observations of the endogenous variable (", 
         length(yVec), ") is not equal to the number of observations", 
         " of the exogenous variables (", nrow(xMat), ")")
  }
  if ("plm.dim" %in% class(data)) {
    dataTable <- matrix(as.integer(data[[1]]), ncol = 1)
    dataTable <- cbind(dataTable, as.integer(data[[2]]))
  }
  else {
    dataTable <- matrix(1:length(yVec), ncol = 1)
    dataTable <- cbind(dataTable, rep(1, nrow(dataTable)))
  }
  nb <- length(xNames)
  dataTable <- cbind(dataTable, yVec)
  if (sum(!is.na(yVec) & is.finite(yVec)) == 0) {
    stop("the dependent variable has no valid observations")
  }
  dataTable <- cbind(dataTable, xMat)
  paramNames <- NULL
  if (nb > 0) {
    for (i in 1:nb) {
      paramNames <- c(paramNames, xNames[i])
      if (sum(!is.na(xMat[, i]) & is.finite(xMat[, i])) == 
          0) {
        stop("regressor '", xNames[i], "' has no valid observations")
      }
    }
  }
  if (is.null(effFormula)) {
    zNames <- NULL
    zIntercept <- FALSE
  }
  else {
    if (class(effFormula) != "formula") {
      stop("argument 'effFormula' must be a formula")
    }
    else if (length(effFormula) != 2) {
      stop("argument 'formula' must be a 1-sided formula")
    }
    me <- match("data", names(mc), 0)
    mfe <- mc[c(1, me)]
    mfe$formula <- effFormula
    attributes(mfe$formula) <- NULL
    mfe$na.action <- na.pass
    mfe[[1]] <- as.name("model.frame")
    mfe <- eval(mfe, parent.frame())
    mte <- attr(mfe, "terms")
    zMat <- model.matrix(mte, mfe)
    if (ncol(zMat) > 0 && colnames(zMat)[1] == "(Intercept)") {
      zIntercept <- TRUE
      zMat <- zMat[, -1, drop = FALSE]
    }
    else {
      zIntercept <- FALSE
    }
    if (nrow(zMat) != nrow(xMat)) {
      stop("the number of observations of the variables explaining", 
           " efficiency (", nrow(zMat), ") is not equal to the number", 
           " of observations of the (regular) regressors (", 
           nrow(xMat), ")")
    }
    dataTable <- cbind(dataTable, zMat)
    zNames <- colnames(zMat)
    if (length(zNames) > 0) {
      for (i in 1:length(zNames)) {
        if (sum(!is.na(zMat[, i]) & is.finite(zMat[, 
                                                   i])) == 0) {
          stop("the regressor for the inefficiency term '", 
               zNames[i], "' has no valid observations")
        }
      }
    }
  }
  nZvars <- length(zNames)
  validObs <- rowSums(is.na(dataTable) | is.infinite(dataTable)) == 
    0
  dataTable <- dataTable[validObs, ]
  nob <- sum(validObs)
  firmId <- sort(unique(dataTable[, 1]))
  nn <- length(firmId)
  firmNo <- rep(NA, nrow(dataTable))
  for (i in 1:nn) {
    firmNo[dataTable[, 1] == firmId[i]] <- i
  }
  dataTable[, 1] <- firmNo
  if (any(is.na(dataTable[, 1]))) {
    stop("internal error: at least one firm number is NA")
  }
  if (min(dataTable[, 1]) != 1) {
    stop("internal error: the smallest firm number must be one")
  }
  if (max(dataTable[, 1]) > nn) {
    stop("internal error: a firm number is larger than the number of firms")
  }
  timeId <- sort(unique(dataTable[, 2]))
  nt <- length(unique(dataTable[, 2]))
  timeNo <- rep(NA, nrow(dataTable))
  for (i in 1:nt) {
    timeNo[dataTable[, 2] == timeId[i]] <- i
  }
  dataTable[, 2] <- timeNo
  if (any(is.na(dataTable[, 2]))) {
    stop("internal error: at least one time period number is NA")
  }
  if (min(dataTable[, 2]) != 1) {
    stop("internal error: the smallest time period number must be one")
  }
  if (max(dataTable[, 2]) > nt) {
    stop("internal error: a time period number is larger", 
         " than the number of time periods")
  }
  for (i in 1:nn) {
    for (j in 1:nt) {
      if (sum(dataTable[, 1] == i & dataTable[, 2] == j) > 
          1) {
        stop("more than one observation for firm '", 
             firmId[i], "' in period '", timeId[j], "'")
      }
    }
  }
  if (modelType == 1) {
    mu <- truncNorm
  }
  else {
    mu <- zIntercept
  }
  if (modelType == 1) {
    eta <- timeEffect
  }
  else {
    eta <- nZvars
  }
  colnames(dataTable) <- c("id", "t", yName, xNames, zNames)
  if (!is.null(rownames(data))) {
    obsNames <- rownames(data)
  }
  else if (!is.null(names(yVec))) {
    obsNames <- names(yVec)
  }
  else if (!is.null(rownames(xMat))) {
    obsNames <- rownames(xMat)
  }
  else if (!is.null(rownames(zMat))) {
    obsNames <- rownames(zMat)
  }
  else {
    obsNames <- NULL
  }
  rownames(dataTable) <- obsNames[validObs]
  names(validObs) <- obsNames
  nParamTotal <- nb + 2 + mu + eta
  if (nParamTotal > nob) {
    stop("the model cannot be estimated,", " because the number of parameters (", 
         nParamTotal, ") is larger than the number of", ifelse(sum(!validObs) > 
                                                                 0, " valid", ""), " observations (", nob, ")")
  }
  if (is.null(startVal)) {
    startVal <- 0
  }
  else {
    if (nParamTotal != length(startVal)) {
      stop("wrong number of starting values (you provided ", 
           length(startVal), " starting values but the model has ", 
           nParamTotal, " parameters)")
    }
  }
  if (nb > 0) {
    ols <- lm(dataTable[, 3] ~ dataTable[, 4:(3 + nb)] - 
                1)
  }
  else if (nb == 0) {
    ols <- lm(dataTable[, 3] ~ -1)
  }
  olsParam <- c(coef(ols), summary(ols)$sigma^2)
  olsStdEr <- sqrt(diag(vcov(ols)))
  olsLogl <- logLik(ols)[1]
  if (nb > 0) {
    gridAdj <- coef(lm(rep(1, nrow(dataTable)) ~ dataTable[, 
                                                           4:(3 + nb)] - 1))
  }
  else {
    gridAdj <- numeric(0)
  }
  if (length(gridAdj) != nb) {
    stop("internal error: the length of 'gridAdj' is not equal to 'nb'.", 
         " Please contact the maintainer of the frontier package")
  }
  returnObj <- .Fortran("front41", modelType = as.integer(modelType), 
                        ineffDecrease = as.integer((!ineffDecrease) + 1), icept = as.integer(0), 
                        nn = as.integer(nn), nt = as.integer(nt), nob = as.integer(nob), 
                        nb = as.integer(nb), mu = as.integer(mu), eta = as.integer(eta), 
                        printIter = as.integer(printIter), indic = as.integer(indic), 
                        tol = as.double(tol), searchTol = as.double(searchTol), 
                        bignum = as.double(bignum), searchStep = as.double(searchStep), 
                        gridDouble = as.integer(gridDouble), gridSize = as.double(gridSize), 
                        maxit = as.integer(maxit), muBound = as.double(muBound), 
                        restartMax = as.integer(restartMax), restartFactor = as.double(restartFactor), 
                        nRestart = as.integer(0), nStartVal = as.integer(length(startVal)), 
                        startVal = as.double(startVal), nRowData = as.integer(nrow(dataTable)), 
                        nColData = as.integer(ncol(dataTable)), dataTable = matrix(as.double(dataTable), 
                                                                                   nrow(dataTable), ncol(dataTable), dimnames = dimnames(dataTable)), 
                        nParamTotal = as.integer(nParamTotal), olsParam = as.double(c(olsParam, 
                                                                                      rep(0, 1 + mu + eta))), gridAdj = as.double(gridAdj), 
                        gridParam = as.double(rep(0, nParamTotal)), startLogl = as.double(0), 
                        mleParam = as.double(rep(0, nParamTotal)), mleCov = matrix(as.double(0), 
                                                                                   nParamTotal, nParamTotal), mleLogl = as.double(0), 
                        nIter = as.integer(0), code = as.integer(0), nFuncEval = as.integer(0))
  if (returnObj$code == 101) {
    stop("the total number of observations exceeds the product of", 
         " the number of firms by the number of years")
  }
  else if (returnObj$code == 102) {
    stop("internal error: calculated variable 'n'", " is not equal to argument 'nParamTotal'.", 
         " Please contact the maintainer of the 'frontier' package", 
         " (arne.henningsen@gmail.com)")
  }
  else if (returnObj$code == 103) {
    stop("wrong number of starting values")
  }
  else if (returnObj$code == 104) {
    stop("a firm number is < 1")
  }
  else if (returnObj$code == 105) {
    stop("a firm number is > number of firms")
  }
  else if (returnObj$code == 106) {
    stop("a period number is < 1")
  }
  else if (returnObj$code == 107) {
    stop("a period number is > number of periods")
  }
  else if (returnObj$code == 108) {
    stop("there are no observations on at least one firm")
  }
  else if (returnObj$code == 109) {
    stop("internal error: 2 + nr - nmu * (im-1)", " is not equal to argument 'nColData'.", 
         " Please contact the maintainer of the 'frontier' package", 
         " (arne.henningsen@gmail.com)")
  }
  else if (returnObj$code > 100) {
    stop("unknown error.", " Please contact the maintainer of the 'frontier' package", 
         " (arne.henningsen@gmail.com)")
  }
  returnObj$nStartVal <- NULL
  returnObj$nRowData <- NULL
  returnObj$nColData <- NULL
  returnObj$nParamTotal <- NULL
  returnObj$ineffDecrease <- as.logical(2 - returnObj$ineffDecrease)
  returnObj$gridDouble <- as.logical(returnObj$gridDouble)
  returnObj$olsParam <- olsParam
  returnObj$olsStdEr <- olsStdEr
  returnObj$olsLogl <- olsLogl
  if (ncol(xMat) == 0) {
    fitVal <- rep(0, sum(validObs))
  }
  else {
    fitVal <- drop(xMat[validObs, , drop = FALSE] %*% returnObj$mleParam[1:nb])
  }
  returnObj$fitted <- matrix(NA, nrow = nn, ncol = nt)
  if (length(fitVal) != nrow(dataTable)) {
    stop("internal error: length of the fitted values is not equal to", 
         " the number of rows of the data table (valid observations)")
  }
  for (i in 1:length(fitVal)) {
    returnObj$fitted[dataTable[i, 1], dataTable[i, 2]] <- fitVal[i]
  }
  resid <- drop(dataTable[, 3] - fitVal)
  returnObj$resid <- matrix(NA, nrow = nn, ncol = nt)
  if (length(resid) != nrow(dataTable)) {
    stop("internal error: length of residuals is not equal to", 
         " the number of rows of the data table (valid observations)")
  }
  for (i in 1:length(resid)) {
    returnObj$resid[dataTable[i, 1], dataTable[i, 2]] <- resid[i]
  }
  returnObj$olsResid <- residuals(ols)
  returnObj$olsSkewness <- moments::skewness(returnObj$olsResid)
  returnObj$olsSkewnessOkay <- returnObj$olsSkewness * (-1)^ineffDecrease >= 
    0
  warnMaxit <- maxit <= returnObj$nIter && maxit > 0
  if (returnObj$mleLogl < returnObj$olsLogl && warnMaxit) {
    warning("the maximum number of iterations has been reached and", 
            " the likelihood value of the ML estimation is less", 
            " than that obtained using OLS;", " please try again using different starting values and/or", 
            " increase the maximum number of iterations")
    warnMaxit <- FALSE
  }
  else if (returnObj$mleLogl < returnObj$olsLogl && maxit > 
           0) {
    warning("the likelihood value of the ML estimation is less", 
            " than that obtained using OLS;", " this indicates that the likelihood maximization did not", 
            " converge to the global maximum or", " that there is no inefficiency", 
            " (you could try again using different starting values)")
  }
  if (warnMaxit) {
    warning("the maximum number of iterations has been reached;", 
            " please try again using different starting values and/or", 
            " increase the maximum number of iterations")
  }
  if (modelType == 1) {
    returnObj$truncNorm <- as.logical(returnObj$mu)
    returnObj$zIntercept <- zIntercept
    returnObj$mu <- NULL
  }
  else {
    returnObj$truncNorm <- truncNorm
    returnObj$zIntercept <- as.logical(returnObj$mu)
    returnObj$mu <- NULL
  }
  if (modelType == 1) {
    returnObj$timeEffect <- as.logical(returnObj$eta)
  }
  else {
    returnObj$timeEffect <- timeEffect
  }
  returnObj$eta <- NULL
  if (returnObj$indic == 2) {
    returnObj$searchScale <- FALSE
  }
  else if (returnObj$indic == 1) {
    returnObj$searchScale <- NA
  }
  else {
    returnObj$searchScale <- TRUE
  }
  returnObj$indic <- NULL
  if (length(startVal) == 1) {
    if (modelType == 1) {
      idx <- 1:(nb + 2)
    }
    else {
      if (nb == 0) {
        idx <- NULL
      }
      else {
        idx <- 1:nb
      }
      idx <- c(idx, (nParamTotal - 1):nParamTotal)
    }
    if (any(returnObj$gridParam[(1:nParamTotal)[-idx]] != 
            0)) {
      warning("internal error: some unused grid-search parameters are", 
              " non-zero: ", paste(returnObj$gridParam, collapse = " "), 
              " please contact the maintainer of the 'frontier' package")
    }
    returnObj$gridParam <- returnObj$gridParam[idx]
    names(returnObj)[names(returnObj) == "startLogl"] <- "gridLogl"
  }
  else {
    returnObj$gridParam <- NULL
  }
  if ("plm.dim" %in% class(data)) {
    rownames(returnObj$resid) <- levels(data[[1]])[firmId]
    colnames(returnObj$resid) <- levels(data[[2]])[timeId]
  }
  else {
    rownames(returnObj$resid) <- obsNames[validObs]
    colnames(returnObj$resid) <- "residuals"
  }
  if (modelType == 2) {
    if (zIntercept) {
      paramNames <- c(paramNames, "Z_(Intercept)")
    }
    if (nZvars > 0) {
      paramNames <- c(paramNames, paste("Z", zNames, sep = "_"))
    }
  }
  if (length(startVal) == 1) {
    returnObj$startVal <- NULL
  }
  paramNames <- c(paramNames, "sigmaSq", "gamma")
  if (modelType == 1) {
    if (truncNorm) {
      paramNames <- c(paramNames, "mu")
    }
    if (timeEffect) {
      paramNames <- c(paramNames, "time")
    }
  }
  if (nb >= 1) {
    betaNames <- paramNames[1:nb]
  }
  else {
    betaNames <- NULL
  }
  names(returnObj$olsParam) <- c(betaNames, "sigmaSq")
  names(returnObj$olsStdEr) <- betaNames
  if (!is.null(returnObj$gridParam)) {
    names(returnObj$gridParam) <- c(betaNames, "sigmaSq", 
                                    "gamma")
  }
  names(returnObj$mleParam) <- paramNames
  rownames(returnObj$mleCov) <- paramNames
  colnames(returnObj$mleCov) <- paramNames
  if (!is.null(returnObj$startVal)) {
    names(returnObj$startVal) <- paramNames
  }
  returnObj$call <- match.call()
  returnObj$validObs <- validObs
  class(returnObj) <- "frontier"
  return(returnObj)
}

estimate_boundary <- function(dat, input, output, sufficient = FALSE,
                              method, AOC = FALSE)
{
  dat$x <- input
  dat$y <- output
  id <- order(dat$x)
  if(method == "Polynomial")
  {
    if(sufficient == TRUE)
    {
      dat$x <- 1 - dat$x
      dat$y <- 1 - dat$y
      odeg <- opt_degree(dat$x, dat$y, dat$x, prange = 0:20)
      boundary <- poly_estimate(dat$x, dat$y, dat$x, deg = odeg)
      boundary < 1 - boundary
    }
    if(sufficient == FALSE)
    {
      odeg <- opt_degree(dat$x, dat$y, dat$x, prange = 0:20)
      boundary <- poly_estimate(dat$x, dat$y, dat$x, deg = odeg)
    }
  }
  if(method == "Kernel")
  {
    if(sufficient == TRUE)
    {
      dat$x <- 1 - dat$x
      dat$y <- 1 - dat$y
      bw <- kern_smooth_bw(dat$x, dat$y)
      boundary <- kernel_smoothing(dat$x, dat$y, dat$x, h = bw)
      boundary <- 1 - boundary
    }
    if(sufficient == FALSE)
    {
      bw <- kern_smooth_bw(dat$x, dat$y)
      boundary <- kernel_smoothing(dat$x, dat$y, dat$x, h = bw)
    }
    
  }
  if(method == "SFA")
  {
    if(sufficient == TRUE)
    {
      dat$x <- 1 - dat$x
      dat$y <- 1 - dat$y
      ka_sfa <- sfa_est(dat$y ~ dat$x | dat$x, 
        ineffDecrease = TRUE, data = dat)
      flip_boundary <- ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * dat$x
      boundary <- 1 - flip_boundary
    }
    if(sufficient == FALSE)
    {
      ka_sfa <- sfa_est(dat$y ~ dat$x | dat$x, 
                                     ineffDecrease = TRUE, data = dat)
      boundary <- ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * input
    }
  }
  if(method == "QR")
  {
    if(sufficient == TRUE)
    {
      dat$x <- 1 - dat$x
      dat$y <- 1 - dat$y
      ka_qr <- rq(dat$y ~ dat$x, tau = 0.95, data = dat)
      flip_boundary <- ka_qr$coefficients[1] + ka_qr$coefficients[2] * dat$x
      boundary <- 1 - flip_boundary
    }
    if(sufficient == FALSE)
    {
      ka_qr <- rq(output ~ input, tau = 0.95, data = dat)
      boundary <- ka_qr$coefficients[1] + ka_qr$coefficients[2] * input
    }
  }
  if(AOC == TRUE)
  {
    id <- order(dat$x)
    divide_by_this <- (max(dat$x) - min(dat$x)) * (max(dat$y) - min(dat$y))
    AOC <- sum(diff(sort(dat$x)) * rollmean(sort(boundary), 2))
    AOC_lower <- sum(diff(dat$x[id]) * min(dat$y))
    AOC_difference <- AOC - AOC_lower
    AOC_percent <- AOC_difference/divide_by_this
    AOC_result <- 1 - AOC_percent
    results_list <- list(boundary, AOC_result)
    names(results_list) <- c("boundary", "AOC")
    return(results_list)
  }
  if(AOC == FALSE)
  {
    return(boundary)
  }
}

AH_AOC <- function(dat, input, output, method, print_style = "ascii",
                   CI = FALSE, nboots)
{
  AOC_estimates <- lapply(method, which_technique_to_use, dat = dat,
                          input = input, output = output)
  estimates <- round(do.call(rbind, AOC_estimates), digits = 3)
  rownames(estimates) <- method
  colnames(estimates) <- "AOC"
  if(CI == TRUE)
  {
    samples <- lapply(method, replicates, nboots = nboots, dat = dat,
                      input = input, output = output)
    CIs <- lapply(samples, quantile, probs = c(0.025, 0.5, 0.975))
    CIs <- round(do.call(rbind, CIs), digits = 3)
    CI_estimates <- cbind(estimates, CIs)
    CI_estimates <- cbind(method, CI_estimates)
    estimates_df <- as.data.frame(CI_estimates)
  }
  if(CI == FALSE)
  {
    add_rows_names <- cbind(method, estimates)
    estimates_df <- as.data.frame(add_rows_names)
  }
  if(print_style == "ascii")
  {
    asciify(estimates_df)
  }
  if(print_style == "R")
  {
    return(estimates_df)
  }
}

which_technique_sim <- function(method, dat, input = input, output = output)
{
  id <- order(input)
  divide_by_this <- (max(input) - min(input)) * (max(output) - min(output))
  ymin <- min(dat$y)
  if(method == "QR"){
    ka_qr <- rq(output ~ input, tau = 0.95, data = dat)
    y_qr <- ka_qr$coefficients[1] + ka_qr$coefficients[2] * input
    bound_function <- function(x = input)
    {
      ka_qr$coefficients[1] + ka_qr$coefficients[2] * x
    }
    integral <- integral(fun = 
      function(x){ka_qr$coefficients[1] + ka_qr$coefficients[2] * x}, 
      xmin = min(input), xmax = max(input))
    qr_result <- (1 - integral)/divide_by_this
    return(qr_result)
  }
  if(method == "SFA"){
    ka_sfa <- sfa_est(output ~ input | input, ineffDecrease = TRUE, data = dat)
    y_sfa <- ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * input
    bound_function <- function(x = input)
    {
      ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * x
    }
    integral <- integral(
      fun = function(x){ka_sfa$mleParam[1] + ka_sfa$mleParam[2] * x}, 
      xmin = min(input), xmax = max(input))
    sfa_result <- (1 - integral)/divide_by_this
    return(sfa_result)
  }
  
  if(method == "Polynomial"){
    odeg <- opt_degree(xtab = input, ytab = output, x = input, prange = 1:20)
    polfront <- poly_estimate(input, output, x = input, deg = odeg)
    AUC_poly <- sum(diff(input[id]) * rollmean(sort(polfront), 2))
    AUC_percent <- AUC_poly/divide_by_this
    poly_result <- (1 - AUC_percent)/divide_by_this
    return(poly_result)
  }
  
  if(method == "Kernel"){
    bw <- kern_smooth_bw(input, output, method = "u", technique = "noh",
                         bw_method="bic")
    kernsmooth <- kernel_smoothing(input, output, input, h = bw)
    AUC_kern <- sum(diff(sort(input)) * rollmean(sort(kernsmooth), 2))
    AUC_percent <- AUC_kern/divide_by_this
    kern_result <- (1 - AUC_kern)/divide_by_this
    return(kern_result)
  }
}

AH_AOC_sim <- function(dat, input, output, method, print_style = "ascii",
                   CI = FALSE, nboots, sim = FALSE)
{
  if(sim == FALSE)
  {
    AOC_estimates <- lapply(method, which_technique_to_use, dat = dat,
                            input = input, output = output)
  }
  if(sim == TRUE)
  {
    AOC_estimates <- lapply(method, which_technique_sim, dat = dat,
                            input = input, output = output)
  }
  estimates <- round(do.call(rbind, AOC_estimates), digits = 3)
  rownames(estimates) <- method
  colnames(estimates) <- "AOC"
  if(CI == TRUE)
  {
    samples <- lapply(method, replicates, nboots = nboots, dat = dat,
                      input = input, output = output)
    CIs <- lapply(samples, quantile, probs = c(0.025, 0.5, 0.975))
    CIs <- round(do.call(rbind, CIs), digits = 3)
    CI_estimates <- cbind(estimates, CIs)
    CI_estimates <- cbind(method, CI_estimates)
    estimates_df <- as.data.frame(CI_estimates)
  }
  if(CI == FALSE)
  {
    add_rows_names <- cbind(method, estimates)
    estimates_df <- as.data.frame(add_rows_names)
  }
  if(print_style == "ascii")
  {
    asciify(estimates_df)
  }
  if(print_style == "R")
  {
    return(estimates_df)
  }
}

